Entanglement at a Two-Dimensional Quantum Critical Point: 
a Numerical Linked Cluster Expansion Study 
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We develop a method to calculate the bipartite entanglement entropy of quantum models, in the 
thermodynamic limit, using a Numerical Linked Cluster Expansion (NLCE) involving only rectan- 
gular clusters. It is based on exact diagonalization of all n x m rectangular clusters at the interface 
between entangled subsystems A and B. We use it to obtain the Renyi entanglement entropy of the 
two-dimensional transverse field Ising model, for arbitrary real Renyi index a. Extrapolating these 
results as a function of the order of the calculation, we obtain universal pieces of the entanglement 
entropy associated with lines and corners at the quantum critical point. They show NLCE to be 
one of the few methods capable of accurately calculating universal properties of arbitrary Renyi 
entropies at higher dimensional critical points. 



Introduction - Quantum critical points (QCPs) [l[ of- 
fer some of the most non-classical, or highly-entangled, 
states in condensed matter physics. Although this high 
degree of entanglement can be a challenge for efforts 
to construct general numerical methods to study QCPs 
@, 9 ; it can also be viewed as a resource to detect and 
classify them. In particular, it is believed that the sub- 
leading scaling terms of the Renyi entanglement entropies 
contain universal coefficients (j, |5[. Thus, these univer- 
sal terms can be studied in quantum many-body mod- 
els using numerical techniques, such as quantum Monte 
Carlo (QMC) or exact diagonalization, and compared 
to quantum field theories as a way of determining the 
universality class of a QCP. In addition to conventional 
universality classes, this procedure has the potential to 
identify unconventional (non-Landau) QCPs, which arc 
predicted to be even more highly entangled than their 
conventional counterparts Thus, it is important to 
have unbiased numerical methods which can calculate 
these universal numbers for a quantitative comparison 
to field theories, models, and perhaps some day, experi- 
mental studies p, @] . 

In this paper, we study the entanglement of a two- 
dimensional (2D) transverse-field Ising model at its quan- 
tum critical point. We use the Numerical Linked Cluster 
Expansion (NLCE) to calculate the Renyi entan- 

glement entropy [1^ for arbitrary index a. The NLCE 
is based on an exact diagonalization of finite-size clus- 
ters up to some "order" O corresponding to the number 
of sites in each cluster. By introducing an innovation 
that allows us to only consider rectangular clusters, we 
are able to perform NLCE up to unprecedented orders. 
Through direct calculation of universal properties associ- 
ated with entanglement across corners, we demonstrate 
that the accuracy of NLCE rivals (or even bests) other 
numerical methods including quantum Monte Carlo, ten- 
sor tree networks, and series expansions. We show that 
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FIG. 1: (a) A sample of cluster shapes and sizes, as used in 
the rectangular NCLE. (b) To calculate the Renyi entropies 
of a given cluster, one considers all divisions into subregions 
A and B defined by translations perpendicular to the line 
(illustrated) or corner. 



the universal term in the entanglement entropy associ- 
ated with a corner is distinct from the value calculated in 
a non-interacting field theory Q . We also use the unique 
ability of NLCE to calculate Renyi entropies with non- 
integer a values to search for a striking change of sign in 
the universal coefficient of line entropy, predicted by an 
interacting field theory at a = 1 [1]. We conclude that 
no sign change takes place, suggesting that either this 
phenomena occurs at exceptionally long length scales, 
or, more hkely, low-order perturbative expansions in the 
field theory are inadequate for the Ising universality class 
in (2 + !)£>. 

Numerical linked cluster expansion (NLCE) [i-Eil IS 



based on expressing an extensive property P of a lat- 
tice model, per site, as a sum over contributions from all 
distinct clusters c that can be embedded in the lattice: 



P/N = L{c) X Wic) 



(1) 
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where L{c) is the number of embeddings of the cluster 
per lattice site, and W{c) is the "weight" of the cluster 
for the property P, defined according to the inclusion- 
exclusion principle, 



W{c) =P(c)-^M^(s). 



(2) 



Here, P{c) is the property calculated for a finite cluster 
c, which contains sub-clusters s. In the NLCE method 
the property is calculated using an exact diagonalization 
for each cluster, and the calculations arc carried out up 
to some order O, the number of sites in the cluster. 

In Eq. ([2]) , the weight of every cluster captures contri- 
butions to P from correlations contained within the size 
of the cluster. When the correlation length grows larger 
than the largest cluster size considered, for example at 
a quantum critical point, results from NCLE require an 
extrapolation in the size of the clusters. Thus, it is im- 
portant to consider as large clusters as possible. 

Typical NLCE approaches involve constructing clus- 
ters that are site- or bond-based. These can involve the 
computationally expensive task of generating all clusters 
c and subclusters s up to a given order O - a proce- 
dure that is related to an NP-complete graph embed- 
ding problem [13|. When calculating groundstatc prop- 
erties, e.g. using Lanczos diagonalization routines, this 
graph embedding problem is the computational bottle- 
neck, restricting the order of conventional expansions to 
^ 16, even with clever use of point-group symmetries and 
topology to reduce graph counting |14| . In this paper, we 
take advantage of the freedom in the NLCE procedure 
to define square-lattice clusters composed of only m x n 
site rectangles (Fig. [Ij, which significantly reduces the 
number of clusters required to self-consistently define the 
sum, and makes the counting of graphs and subgraphs 
trivial [l5|, 16 1. The computational bottleneck becomes 
the Lanczos diagonalization of each rectangular cluster, 
allowing us to push the NLCE to significantly higher or- 
ders - up to O = 26 with moderate effort on simple 
desktop workstations [i^. We stress that the NLCE is 
not a calculation for a finite size system, rather a system- 
atic approximation for the thermodynamic limit, where 
the rectangles provide a way to sum up contributions 
from different length scales corresponding to the order of 
the cluster. NLCE can systematically encapsulate signif- 
icantly larger-range correlations than conventional finite- 
size studies of toroidal clusters (as is often done in Lanc- 
zos diagonalizations) which are crucial to determine the 
singular behavior at a QCP [iyj . 

We use this method to study the scaling of the gen- 
eralized Renyi entropies at the quantum critical point of 
the 2D transverse field Ising model (TFIM), 



H 



(3) 



where '^i is a Pauli spin operator, so that cr^ has eigen- 
values ±1. In this equation, the first sum is over lat- 
tice bonds, while the second sum is over lattice sites. 
To calculate the Renyi entanglement entropies, Sa = 
ln(Trp^)/(I — a) one must position each cluster in rela- 
tion to the boundary between subregions A and B, and 
calculate the reduced density matrix for the subsystem 
A. Translational symmetry along the boundary simply 
gives an overall factor of length L, and automatically pro- 
duces the entropy per unit length associated with a line, 
when only translationally distinct clusters are included. 
To take care of translations perpendicular to the line, or 
translations with respect to a corner, it is equivalent and 
more convenient to consider a given n x m cluster only 
once, but allow all possible linear (in the case of a line, 
illustrated in Fig. [T{b)) or quadrant-based divisions of 
that cluster. 

A key advantage of the NLCE method over quantum 
Monte Carlo (QMC) [l3| and series expansions [31 is 
that it can be used to calculate Renyi entropies for any a 
value including a < 1. Lanczos diagonalization provides 
us with the exact ground state wavefunction of each clus- 
ter, allowing an explicit calculation of the reduced den- 
sity matrix, from which Renyi entropies with arbitrary 
a can be obtained. Another distinct advantage of the 
NLCE method over QMC is that one can analytically 
separate the Renyi entropies associated with lines and 
corners. When the subsystem ^ is a half-plane, one ob- 
tains only the entropy associated with the line. When the 
subsystem A is a quadrant, it contains both line and cor- 
ner contributions. A suitable choice of subdivision of the 
system into half-planes and corners is sufHcient to isolate 
the corner contribution from every graph, thus leading 
to a separate calculation for line and corner entropies for 
the TFIM, that we call Sa = Sa/L and Cq, below. This, in 
turn, allows a more accurate determination of the singu- 
larities associated with each term than possible in QMC, 
where e.g. the dominant "area law" can easily overwhelm 
sub-leading terms such as corner contributions. 

In Fig. [2] we show results for line and corner entropies 
for the TFIM as a function of h/J at different orders, 
where hc/J~ 3.044 is the quantum critical point. Prcvi- 
ousW, the series expansion method was used to calculate 
S2 [18[ , but not for example 5*1 , which is also inaccessible 
to QMC calculations due to their reliance on the replica 
trick. Note that, in order to get convergence for h < he, 
one needs to add a static ordered moment for the sites 
outside the cluster. These moments apply a boundary 
field on the sites in the cluster, via the exchange J, and 
ensure that fluctuations remain bounded. This is analo- 
gous to low field series expansions, and we call it low field 
NLCE. For h > he no such boundary field is needed and 
we call it high field NLCE. 

The rest of the paper will focus on the quantum critical 
behavior. Since high field NLCE is significantly more 
accurate, we will restrict our study to h > he- General 
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FIG. 2: The corner term C2 and line term S2 as a function 
of h/J for the second Renyi entanglement entropy using both 
high field and low field NLCE of orders 8, 12, 16, and 24. The 
dotted line denotes he/ J = 3.044 and the dashed line shows 
the series expansion results [l8| . 



arguments and our numerical study show that one can 
associate a length scale with order O which goes as VO. 
In Fig. 131 we show the behavior of the corner term Cq 
at h = he- It is predicted to scale like Ca ln(L) 
[3], with Qa being universal. As shown in the inset of 
Fig.m a plot of Ca versus ln(l/VO) can be extrapolated 
to the O — > oo limit to obtain the universal term Oq. 
In the main plot of Fig. [31 we show Oq, as a function of 
Renyi index. The values of this universal constant, as 
calculated by other methods, are also shown on the plot. 

The coefficient 02 has been calculated for the 2D TFIM 
using several different numerical methods. In Ref. [l8| . 
series expansion found 02 = —0.0055(5) [l^. Finite- 
temperature QMC calculations on a single 36 x 36-size 
lattice report -0.0075(25) Our current NLCE calcu- 
lation gives —0.0053, with an uncertainty in the last digit 
- a number consistent with the series expansion result. 

This value of 02 for the TFIM has been compared sev- 
eral times in the past literature to the value calculated 
analytically by Casini and Huerta for a free scalar field 
theory, 02 = -0.0064 [llj. We note that, ahhough this 
is numerically close to our value of —0.0053, in fact one 
should not expect correspondence since this is an inter- 
acting model that is not described by a free Gaussian 
fixed point. Rather, one needs a calculation of the quan- 
tity in an interacting theory - the lack of which we hope 
will motivate future field-theoretic calculations of da at 
the Wilson-Fisher fixed point. 

To our knowledge, the universal term ai has only been 
calculated once previously using a tensor network vari- 
ational ansatz called the "tensor tree network" (TTN), 
by Tagliacozzo and co-workers Their value for the 
universal coefficient is oi = —0.0095(1). The value from 
free scalar field theory is —0.012 The value from the 
present NLCE work up to order O = 26 is ai = —0.0140, 
a number that is relatively close to the free field value. 



FIG. 3: The coefficient —a^ from fitting to the function 
aa log y/O + ba (purple dots) plotted along with other numer- 
ical estimates of —aa from free scalar field theory 01 (green 
diamonds), tensor tree network [22| (red star), finite-T QMC 
[20I I (yellow square), and series expansion O, (3 (black cir- 
cle). Inset: Ca for a = 1, 2, 10 along with fits plotted vs 1/VO 
on a logarithmic scale. 



Finally we note, based on the shape of the aa versus 
a curve in relation to the results from non-interacting 
theory in Fig. |3l it is possible that a — 2 may not be 
the optimal Renyi index to make comparisons. This is 
particularly important for methods such as QMC or se- 
ries expansion that are restricted to integer a > 2 due to 
reliance on the replica trick. In those cases, larger alpha 
values (e.g. 5*3) might be more suitable to distinguish the 
free scalar field theory from interacting ones. 

Next, we examine the behavior of the line term s„. 
Near the quantum critical point the line entropy is ex- 
pected to take the form [5|] , 



S2=vL + SoiL/0, 



(4) 



where, the coefficient rj is non-universal but Sq is a uni- 
versal function of its argument. In the limit L/S. — >■ 0, it 
becomes a universal number 7 associated with the criti- 
cal point, while for L/^ — > 00, it takes the form r^L/^, 
with a universal amplitude ratio Tq,. Like Ca, this num- 
ber may be used to distinguish diflFerent critical points 
by comparison between field theory, and numerical cal- 
culation on model systems. For the Gaussian fixed point, 
TqCx- (1 + ^) 0, [i^l- More interestingly, for the inter- 
acting field theory relevant to the TFIM, Ref. Q reports 
ra oc a — 1, a result consistent between e and 1/N ex- 
pansions at the Wilson- Fisher fixed point. According to 
this, the universal amplitude ratio should change sign as 
a function of Renyi index a, as one passes through the 
von Neumann point a ^ 1. 

To obtain ra, we first need to obtain Sa in the ther- 
modynamic limit for h near he. We expect partial sums 
to order O to behave as 1/VO if the correlation length 
of the system ^ ^ \h — hc\~'^ is larger than \/0, but to 
saturate to the thermodynamic value when ^ is less than 
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\fO . Hence, the NLCE data for h near are fit for each 
h and a to the function, 



0.3 



(5) 
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The constant C was obtained from a 
Even if C is aUowed 



with 

fit to series expansion data [18 
to vary in the NLCE fits, it does not change significantly 
as a function of a. Thus it was held fixed at this value for 
all fits to Eq. ^ , a few of which are shown in Fig. IH^a) 
for si at three different values of h. 

Figure HJ^b) shows an example of the resulting linear 
term, A as a function oi h ioi a — 1, from the fits to 
Eq. (O. This term is then fit to, 



A{h) = D + ro,\h - hcl' 



(6) 



for all values of a. This gives Fig. IH^c), showing the 
ratio fjj = Ta/ri obtained after eliminating an unknown 
a independent prefactor in Tq. It is apparent from this 
figure that there is no singularity or sign change in the 
universal amplitude Tq at the von Neumann point a = 1. 

To independently confirm this result, we have used an 
alternative simpler fitting procedure to obtain Tq, using 
data only at the critical point h^. As discussed above, 
the quantity y/O is a measure of the correlation length ^ 
explored (see Fig. 1) by NLCE. Hence, one can conjecture 
that the slope in the plot of Sa versus \/0 is proportional 
to Va at he- Since there is an unknown proportionality 
constant between y/O and ^, one can examine the ratio 
of the slope for a given Renyi index a with the slope 
at a = 1 to eliminate that constant. The result of this 
procedure gives similar curves as in Fig.|4l^c) for a > 0.5, 
with a slightly stronger divergence of Tq, for a < 0.5. 
Again, although unusual non-monotonic behavior occurs 
around a ~ 1, there is no singularity or sign change of 
Tq, at the von Neumann point. 

Discussion - In this paper we have developed a new 
NLCE procedure for calculating the bipartite Renyi en- 
tanglement entropy of quantum lattice models using all 
nx m rectangular clusters at the interface of subsystems 
A and B. The use of rectangular clusters, which can 
straight-forwardly be modified for other two-dimensional 
lattices, overcomes the computational bottleneck aris- 
ing from the exponentially-scaling subgraph isomorphism 
problem that restricted previous NLCE methods to order 
O < 16. With this innovation, our NLCE method is only 
constrained by the exact diagonalization kernel - easily 
obtaining results up to order 26 in the current study. 

Using the NLCE method, we have calculated individ- 
ually the line and corner contributions of the Renyi en- 
tanglement entropy of the transverse field Ising model 
for arbitrary Renyi index a. Extrapolating in the order 
of the calculation allows us to estimate universal critical 
properties of this model. For entanglement across cor- 
ners, we conclusively demonstrate that a universal term 
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FIG. 4: (a) Fits of the NLCE line term s to Eq. Q for a = 1 
from NLCE for h = 3.048(top), 3.075, and 3.1(bottom). (b) 
The resulting linear term A{h) from the fits in inset (a), fit 
to Eq. (|6]). (c) The normalized coefficient fa as a function of 



in the a = 2 entropy is distinct from the value calculated 
in a non-interacting field theory by Casini and Huerta 
We hope this motivates future field-theory calcula- 
tions at the interacting fixed point. For entanglement 
across lines, we searched for the striking sign change in 
the universal coefficient predicted by an interacting field 
theory at a = 1 0], but find that no such sign change 
takes place, possibly implying a need for higher order e 
or 1 /N expansion terms in the field theory. 

We have demonstrated that the accuracy of NLCE for 
the extraction of universal critical entanglement prop- 
erties rivals other numerical methods including quantum 
Monte Carlo, tensor tree networks, and series expansions. 
In addition to the simplicity in graph counting introduced 
here, this success also stems from the fact that NLCE 
can calculate Renyi entanglement entropies for arbitrary 
index a, and analytically separate contributions to the 
entanglement from line and corner geometries. 

Perhaps the most appealing feature of the method is 
that this accuracy results from a very simple numerical 
algorithm - less than a thousand lines of code are required 
to implement the rectangular-cluster NLCE into any ex- 
isting Lanczos program. Rectangular cluster geometries 
also give the exciting option of using the density ma- 
trix renormalization group (25j as a cluster solver, which 
could increase the accessible order of the NLCE to 100 or 
beyond. Ultimately, its high accuracy, coupled with sim- 
plicity of implementation, lends hope that NLCE will be 
widely adopted to the study of entanglement and other 
universal properties at quantum critical points in two and 
higher dimensions in the future. 
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